clc; clear all;
format long
%% 生成仿真信号波形，sig_nonoise为无噪声信号(原信号--目标信号)，sig_noise为添加噪声后的信号（加噪信号）
rng(123); %设置随机种子，保证每次生成噪声的一致性

% %生成高频信号
% Fs=100; %数据采样率Hz
% t=(1:1/Fs:2048*1/Fs)'; %对数据进行采样,将t转置为1列--生成时间域
% N = length(t); %设置时间域范围：0-N（连续信号）N：数据的采样数目
% f1 =0.8; %信号的频率
% f2=0.05;
% sig_nonoise=2*sin(2*pi*f1*t+cos(2*pi*f2*t)); %生成原始信号
% % sig_nonoise = sig_nonoise - sign(t-0.3) - sign(0.72-t);%生成原信号
% sig_noise = sig_nonoise + 0.2*randn(size(t));%生成加噪信号

%生成低频信号
% N = 201;
% t = linspace(0,1,N);
% sig_nonoise = 4*sin(4*pi*t);
% sig_nonoise = sig_nonoise - sign(t-0.3) - sign(0.72-t);
% sig_noise = sig_nonoise + 0.2*randn(size(t));

% plot(t,sig)
% 生成添加高斯白噪声的加噪信号
% nt=0.9*randn(N,1); %高斯白噪声生成
% sig_noise = sig_nonoise + nt;%生成加噪信号

%% 模拟仿真：录制音频处理
R = audiorecorder(44100,16,1);%参数：采样频率，存储位数，通道数；
%录音
% record(R);
% stop(R);
% myword=getaudiodata(R);
% plot(myword);
% audiowrite('D:\Mytool\Matlab 2020a\bin\souds.wav',myword,44100);
%  [y,fs]=audioread('D:\Mytool\Matlab 2020a\bin\souds.wav'); 
[y,fs]=audioread('H:\测试用\博学楼59号楼 2 副本.wav'); %用自己存放的测试用音频地址
sound(y,fs);
y = y(107501:117501);
n=length(y);
y_f=fft(y,n);
f=fs*(1:n/2)/n;
subplot 211;
plot(y);
xlabel('时间s');
ylabel('幅值 ');
title('加噪前的时域波形');
subplot 212;
plot(f,abs(y_f(1:n/2,1)));
xlabel('频率Hz');
ylabel('频率幅值');
title('加噪前的频谱图');
%对原始的语音信号进行加噪处理，
% 使用randn函数产生一个与音频信号等长度的高斯随机噪声信号（噪声年夜小取决于随机函数的幅度倍数）。
% 然后通过信号的叠加y_z=y+noise;，生成一个新的信号。再通过傅里叶变换获得信号在频域上的波形，
% 最后画出时域与频域波形。
L=length(y);
noise=0.1*randn(L,1);
y_z=y+noise;
% sound(y_z,fs);
n=length(y);
y_zf=fft(y_z,n);
f=fs*(1:n/2)/n;
subplot 211;
plot(y_z);
xlabel('时间s');
ylabel('幅值');
title('加噪后的时域波形');
subplot 212;
plot(f,abs(y_zf(1:n/2)));
xlabel('频率Hz');
ylabel('频率幅值');
title('加噪后的频谱图');
sig_nonoise = y;
sig_noise = y_z;
N=length(sig_nonoise);

%% 初始参数的记录
%iceemdan参数
flag = 1;% 当flag为1时获得理想信号，当flag为0时循环进行iceemdan
expect = 5;% 期望，预期误差标准差的值，越结果越精确，但不宜过小（会导致迭代次数过多或无法得出结果）一般情况均大于1，小于10
Nstd = 0.2;% Nstd为附加噪声标准差与Y标准差之比
NR = 1;
NE = 100;% NE为对信号的平均次数
MaxIter = 5000;% MaxIter 最大迭代次数
expect_p = 0.2;%取值越大，精度越高（取值过大会消除重要特征）%高频仿真：0.2，低频仿真：0.4
expect_per = 0.002;%取值越小，经度越高（取值过小会消除重要特征）%高频仿真：0.002，低频仿真：0.08
%小波阈值去噪参数
Psig=sum(sig_nonoise'*sig_nonoise)/N;
Pnoi1=sum((sig_noise-sig_nonoise)'*(sig_noise-sig_nonoise))/N;
SNR1=10*log10(Psig/Pnoi1);%生成的初始加噪信号的信噪比
%构建结构体option
% option = struct('N',{},'a4_alpha',{},'a4_gamma',{},'a4_p',{},'a5_alpha',{},'a5_beta',{},'a6_alpha',{},'a6_beta',{});
option.N = N;

%% 调用iceemdan进行去噪
%% ICEEMDAN
ecg = sig_noise';%用于高频信号
% ecg = sig_noise;%用于低频信号
% ecg = sig';
while flag
    [modes]=iceemdan(ecg,Nstd,NR,MaxIter,1);%iceemdan
    modes=modes';
    t=1:length(ecg);
    [a,b]=size(modes);
    for i=2:a
        if std(modes(i-1,:)-sig_nonoise) < expect
            flag = 0;
        end
    end
end
per = getper(modes,ecg,a);
modes_all = 0;
for i = 2:1:a-1
    [r,p] = corr(modes(i,:)',ecg','type','Pearson');
    p_target(i) = p;
%     disp(per(i));
%     disp(p_target(i));
    if per(i)<expect_per || p_target(i)>expect_p
        modes_all = modes_all + modes(i,:);
    end
    Result = modes_all;
end
% sig_noise = Result';%赋值转换
% sig = Result';%赋值转换

%% 调用函数进行滤波
wave_name = 'db4';  %小波名称
wave_level = 3;        %小波分解水平
threshold = 's6';     %阈值类型
% SORH_1 = 's4';     %阈值1类型
% SORH_2 = 's5';     %阈值2类型
%a4阈值函数的参数
option.a4_alpha = 2;
option.a4_gamma = 0.5;
option.a4_p = 0.2;%p=0-thr
%a5阈值函数的参数
option.a5_alpha = 1;
option.a5_beta = 0.6;
%a6阈值函数的参数
option.a6_alpha = 0.5;
option.a6_beta = 0.5;
option.a6_sigma = 1;
option.a6_mu = 0;
option.a6_gamma = 0.0405;%取值最好在0.035-0.045（最佳值在0.0405附近）
option.a6_per = 0.7;%小波阈值去噪的权重（调节权重来调节最终复原波的小波阈值去噪与ICEEMDAN的占比，以此适应不同情况）
% select_topic = 'sqtwolog'; %阈值选择规则
select_topic = 'rigrsure'; %阈值选择规则
%调用复合阈值函数对信号进行降噪处理
% sig1 = filterWaveletTh(sig,wname,SORH_1,lev,tptr,option);  %调用函数进行滤波
% s = filterWaveletTh(sig1,wname,SORH_2,lev,tptr,option);  %调用函数进行滤波
%调用单个阈值函数对信号进行降噪处理
if strcmp(threshold ,'h')
    sig = HardThreshold(sig_noise,wave_name,wave_level,select_topic,option);  %调用函数进行滤波
    disp('启用硬阈值函数');
elseif strcmp(threshold ,'s')
    sig = SoftThreshold(sig_noise,wave_name,wave_level,select_topic,option);  %调用函数进行滤波
    disp('启用软阈值函数');
elseif strcmp(threshold ,'s1')
    sig = S1Threshold(sig_noise,wave_name,wave_level,select_topic,option);  %调用函数进行滤波
    disp('启用s1阈值函数');
elseif strcmp(threshold ,'s2')
    sig = S2Threshold(sig_noise,wave_name,wave_level,select_topic,option);  %调用函数进行滤波
    disp('启用s2阈值函数');
elseif strcmp(threshold ,'s3')
    sig = S3Threshold(sig_noise,wave_name,wave_level,select_topic,option);  %调用函数进行滤波
    disp('启用s3阈值函数');
elseif strcmp(threshold ,'s4')
    sig = S4Threshold(sig_noise,wave_name,wave_level,option);  %调用函数进行滤波
    disp('启用s4阈值函数');
elseif strcmp(threshold ,'s5')
    sig = S5Threshold(sig_noise,wave_name,wave_level,select_topic,option);  %调用函数进行滤波
    disp('启用s5阈值函数');
elseif strcmp(threshold ,'s6')
    sig = S6Threshold(sig_noise,wave_name,wave_level,option);  %调用函数进行滤波
    disp('启用s6阈值函数');
end

%% 综合双去噪
sig_result = option.a6_per*sig+(1-option.a6_per)*Result';

%% 对于去噪效果的评估指标
sig_origin = sig_nonoise;  %无噪声信号
sig_filter = sig_result;  %滤波后信号
[SNR,MSE,NCC] = FilterEffectEvaluation(sig_origin,sig_filter,N);
% 对滤波效果进行评估
% 目前包含的评估指标有：信噪比SNR，均方差MSE，波形相似系数NCC
% 输入：
% sig_origin：无噪声的原始数据，一维序列
% sig_filter：滤波后的数据，一维序列
disp(['滤波后的SNR、MSE和NCC值分别为：',num2str(SNR),'、',num2str(MSE),'、',num2str(NCC)]);

%% 画滤波前后对比图
figure('color','w')
subplot(311);plot(sig_nonoise,'k');title('原始信号（未加入噪声）')
subplot(312);plot(sig_noise,'k');title('原始信号（加入噪声）')
subplot(313);plot(sig,'k');title('滤波后信号')
% 注意，对于不知道纯净信号的情况，无法绘制“未加入噪声”的信号图像
% 此时需要将上述4行代码删除，并使用下述3行代码：
% figure('color','w')
% subplot(211);plot(sig_noise,'k');title('原始信号')
% subplot(212);plot(sig,'k');title('滤波后信号')

%% 2.ICEEMDAN分解图
figure;
subplot(a,1,1);
plot(t,ecg);% the ECG signal is in the first row of the subplot
ylabel('original');
set(gca,'xtick',[]);
title('ICEEMDAN');
axis tight;
for i=2:a
    subplot(a+1,1,i);
%     disp(std(modes(i-1,:)-sig))
    plot(t,modes(i-1,:));
    ylabel (['IMF ' num2str(i-1)]);
    set(gca,'xtick',[]);
    xlim([1 length(ecg)]);
end

subplot(a+1,1,a+1);
plot(t,Result);
ylabel('RES ');
xlim([1 length(ecg)]);
xlabel('样本点');
% 画信号ICEEMDAN分解图
% 输入：
% data为待分解信号
% FsOrT为采样频率或采样时间向量，如果为采样频率，该变量输入单个值；如果为时间向量，该变量为与y相同长度的一维向量。如果未知采样频率，可设置为1
% Nstd为附加噪声标准差与Y标准差之比
% NE为对信号的平均次数
% MaxIter：最大迭代次数
% 输出：
% imf为经ICEEMDAN分解后的各imf分量值

%% 实现小波阈值滤波的filterWaveletTh
%% 阈值函数1 -- 硬阈值 -- HardThreshold
function sig = HardThreshold(data,wave_name,wave_level,select_topic,option)

    %% 用wave_name小波对含噪信号进行wave_level层分解并提取系数
    [c,l]=wavedec(data,wave_level,wave_name);
    %取第wave_level层低频近似系数
    ca=appcoef(c,l,wave_name,wave_level);
    %取各层高频细节系数
    thr=thselect(data,select_topic); % 阈值获取
    if wave_level==1
        cd1=detcoef(c,l,1);
        %% 进行threshold阈值处理
        ysoft1=wthresh(cd1,'h',thr);
        c1 = [ca;ysoft1];
    end
    if wave_level==2
        cd1=detcoef(c,l,1);
        cd2=detcoef(c,l,2);
        %% 进行threshold阈值处理
        ysoft1=wthresh(cd1,'h',thr);
        ysoft2=wthresh(cd2,'h',thr);
        c1 = [ca;ysoft2;ysoft1];
    end
    if wave_level==3
        cd1=detcoef(c,l,1);
        cd2=detcoef(c,l,2);
        cd3=detcoef(c,l,3);
        %% 进行threshold阈值处理
        ysoft1=wthresh(cd1,'h',thr);
        ysoft2=wthresh(cd2,'h',thr);
        ysoft3=wthresh(cd3,'h',thr);
        c1 = [ca;ysoft3;ysoft2;ysoft1];
    end
    if wave_level==4
        cd1=detcoef(c,l,1);
        cd2=detcoef(c,l,2);
        cd3=detcoef(c,l,3);
        cd4=detcoef(c,l,4);
        %% 进行threshold阈值处理
        ysoft1=wthresh(cd1,'h',thr);
        ysoft2=wthresh(cd2,'h',thr);
        ysoft3=wthresh(cd3,'h',thr);
        ysoft4=wthresh(cd4,'h',thr);
        c1 = [ca;ysoft4;ysoft3;ysoft2;ysoft1];
    end
    if wave_level==5
        cd1=detcoef(c,l,1);
        cd2=detcoef(c,l,2);
        cd3=detcoef(c,l,3);
        cd4=detcoef(c,l,4);
        cd5=detcoef(c,l,5);
        %% 进行threshold阈值处理
        ysoft1=wthresh(cd1,'h',thr);
        ysoft2=wthresh(cd2,'h',thr);
        ysoft3=wthresh(cd3,'h',thr);
        ysoft4=wthresh(cd4,'h',thr);
        ysoft5=wthresh(cd5,'h',thr);
        c1 = [ca;ysoft5;ysoft4;ysoft3;ysoft2;ysoft1];
    end
    sig=waverec(c1,l,wave_name);
    disp(thr);
end
%% 阈值函数2 -- 软阈值 -- SoftThreshold
function sig = SoftThreshold(data,wave_name,wave_level,select_topic,option)
     %% 用wave_name小波对含噪信号进行wave_level层分解并提取系数
    [c,l]=wavedec(data,wave_level,wave_name);
    %取第wave_level层低频近似系数
    ca=appcoef(c,l,wave_name,wave_level);
    %取各层高频细节系数
    thr=thselect(data,select_topic); % 阈值获取
%调用matlab自带函数    
    if wave_level==1
        cd1=detcoef(c,l,1);
        %% 进行threshold阈值处理
        ysoft1=wthresh(cd1,'s',thr);
        c1 = [ca;ysoft1];
    end
    if wave_level==2
        cd1=detcoef(c,l,1);
        cd2=detcoef(c,l,2);
        %% 进行threshold阈值处理
        ysoft1=wthresh(cd1,'s',thr);
        ysoft2=wthresh(cd2,'s',thr);
        c1 = [ca;ysoft2;ysoft1];
    end
    if wave_level==3
        cd1=detcoef(c,l,1);
        cd2=detcoef(c,l,2);
        cd3=detcoef(c,l,3);
        %% 进行threshold阈值处理
        ysoft1=wthresh(cd1,'s',thr);
        ysoft2=wthresh(cd2,'s',thr);
        ysoft3=wthresh(cd3,'s',thr);
%         c1 = [ca;ysoft3';ysoft2';ysoft1'];%用于高频信号
        c1 = [ca;ysoft3;ysoft2;ysoft1];%用于高频信号
%         c1 = [ca,ysoft3,ysoft2,ysoft1];%用于低频信号
    end
    if wave_level==4
        cd1=detcoef(c,l,1);
        cd2=detcoef(c,l,2);
        cd3=detcoef(c,l,3);
        cd4=detcoef(c,l,4);
        %% 进行threshold阈值处理
        ysoft1=wthresh(cd1,'s',thr);
        ysoft2=wthresh(cd2,'s',thr);
        ysoft3=wthresh(cd3,'s',thr);
        ysoft4=wthresh(cd4,'s',thr);
        c1 = [ca;ysoft4;ysoft3;ysoft2;ysoft1];
    end
    if wave_level==5
        cd1=detcoef(c,l,1);
        cd2=detcoef(c,l,2);
        cd3=detcoef(c,l,3);
        cd4=detcoef(c,l,4);
        cd5=detcoef(c,l,5);
        %% 进行threshold阈值处理
        ysoft1=wthresh(cd1,'s',thr);
        ysoft2=wthresh(cd2,'s',thr);
        ysoft3=wthresh(cd3,'s',thr);
        ysoft4=wthresh(cd4,'s',thr);
        ysoft5=wthresh(cd5,'s',thr);
        c1 = [ca;ysoft5;ysoft4;ysoft3;ysoft2;ysoft1];
    end
    sig=waverec(c1,l,wave_name);
    disp(thr);

%调用自定义函数
%     if wave_level==1
%         cd1=detcoef(c,l,1);
%         %% 进行threshold阈值处理
%         ysoft1=SGetThreshold(cd1,option,thr);
%         c1 = [ca;ysoft1'];
%     end
%     if wave_level==2
%         cd1=detcoef(c,l,1);
%         cd2=detcoef(c,l,2);
%         %% 进行threshold阈值处理
%         ysoft1=SGetThreshold(cd1,option,thr);
%         ysoft2=SGetThreshold(cd2,option,thr);
%         c1 = [ca;ysoft2';ysoft1'];
%     end
%     if wave_level==3
%         cd1=detcoef(c,l,1);
%         cd2=detcoef(c,l,2);
%         cd3=detcoef(c,l,3);
%         %% 进行threshold阈值处理
%         ysoft1=SGetThreshold(cd1,option,thr);
%         ysoft2=SGetThreshold(cd2,option,thr);
%         ysoft3=SGetThreshold(cd3,option,thr);
%         c1 = [ca;ysoft3';ysoft2';ysoft1'];
%     end
%     if wave_level==4
%         cd1=detcoef(c,l,1);
%         cd2=detcoef(c,l,2);
%         cd3=detcoef(c,l,3);
%         cd4=detcoef(c,l,4);
%         %% 进行threshold阈值处理
%         ysoft1=SGetThreshold(cd1,option,thr);
%         ysoft2=SGetThreshold(cd2,option,thr);
%         ysoft3=SGetThreshold(cd3,option,thr);
%         ysoft4=SGetThreshold(cd4,option,thr);
%         c1 = [ca;ysoft4';ysoft3';ysoft2';ysoft1'];
%     end
%     if wave_level==5
%         cd1=detcoef(c,l,1);
%         cd2=detcoef(c,l,2);
%         cd3=detcoef(c,l,3);
%         cd4=detcoef(c,l,4);
%         cd5=detcoef(c,l,5);
%         %% 进行threshold阈值处理
%         ysoft1=SGetThreshold(cd1,option,thr);
%         ysoft2=SGetThreshold(cd2,option,thr);
%         ysoft3=SGetThreshold(cd3,option,thr);
%         ysoft4=SGetThreshold(cd4,option,thr);
%         ysoft5=SGetThreshold(cd5,option,thr);
%         c1 = [ca;ysoft5';ysoft4';ysoft3';ysoft2';ysoft1'];
%     end
%     sig=waverec(c1,l,wave_name);
%     disp(thr);
end
    
%% 阈值函数3 -- 改进方法s1 -- S1Threshold  
function sig = S1Threshold(data,wave_name,wave_level,select_topic,option)
     %% 用wave_name小波对含噪信号进行wave_level层分解并提取系数
     [c,l]=wavedec(data,wave_level,wave_name);
     %disp(c);
     %取第wave_level层低频近似系数
     ca=appcoef(c,l,wave_name,wave_level);
     %取各层高频细节系数
     thr=thselect(data,select_topic); % 阈值获取
     if wave_level==1
        cd1=detcoef(c,l,1);
        %% 进行threshold阈值处理
        ysoft1=S1GetThreshold(cd1,option,thr);
        c1 = [ca;ysoft1'];
    end
    if wave_level==2
        cd1=detcoef(c,l,1);
        cd2=detcoef(c,l,2);
        %% 进行threshold阈值处理
        ysoft1=S1GetThreshold(cd1,option,thr);
        ysoft2=S1GetThreshold(cd2,option,thr);
        c1 = [ca;ysoft2';ysoft1'];
    end
    if wave_level==3
        cd1=detcoef(c,l,1);
        cd2=detcoef(c,l,2);
        cd3=detcoef(c,l,3);
        if cd3 > thr
        disp(cd3);
        end
        %% 进行threshold阈值处理
        ysoft1=S1GetThreshold(cd1,option,thr);
        ysoft2=S1GetThreshold(cd2,option,thr);
        ysoft3=S1GetThreshold(cd3,option,thr);
        c1 = [ca;ysoft3';ysoft2';ysoft1'];%用于高频信号
%         c1 = [ca,ysoft3,ysoft2,ysoft1];%用于低频信号
    end
    if wave_level==4
        cd1=detcoef(c,l,1);
        cd2=detcoef(c,l,2);
        cd3=detcoef(c,l,3);
        cd4=detcoef(c,l,4);
        %% 进行threshold阈值处理
        ysoft1=S1GetThreshold(cd1,option,thr);
        ysoft2=S1GetThreshold(cd2,option,thr);
        ysoft3=S1GetThreshold(cd3,option,thr);
        ysoft4=S1GetThreshold(cd4,option,thr);
        c1 = [ca;ysoft4';ysoft3';ysoft2';ysoft1'];
    end
    if wave_level==5
        cd1=detcoef(c,l,1);
        cd2=detcoef(c,l,2);
        cd3=detcoef(c,l,3);
        cd4=detcoef(c,l,4);
        cd5=detcoef(c,l,5);
        %% 进行threshold阈值处理
        ysoft1=S1GetThreshold(cd1,option,thr);
        ysoft2=S1GetThreshold(cd2,option,thr);
        ysoft3=S1GetThreshold(cd3,option,thr);
        ysoft4=S1GetThreshold(cd4,option,thr);
        ysoft5=S1GetThreshold(cd5,option,thr);
        c1 = [ca;ysoft5';ysoft4';ysoft3';ysoft2';ysoft1'];
    end
    sig=waverec(c1,l,wave_name);
    disp(thr);
end

%% 阈值函数4 -- 改进方法s2 -- S2Threshold  
function sig = S2Threshold(data,wave_name,wave_level,select_topic,option)
     %% 用wave_name小波对含噪信号进行wave_level层分解并提取系数
     [c,l]=wavedec(data,wave_level,wave_name);
     %disp(c);
     %取第wave_level层低频近似系数
     ca=appcoef(c,l,wave_name,wave_level);
     %取各层高频细节系数
     thr=thselect(data,select_topic); % 阈值获取
     if wave_level==1
        cd1=detcoef(c,l,1);
        %% 进行threshold阈值处理
        ysoft1=S2GetThreshold(cd1,option,thr);
        c1 = [ca;ysoft1'];
    end
    if wave_level==2
        cd1=detcoef(c,l,1);
        cd2=detcoef(c,l,2);
        %% 进行threshold阈值处理
        ysoft1=S2GetThreshold(cd1,option,thr);
        ysoft2=S2GetThreshold(cd2,option,thr);
        c1 = [ca;ysoft2';ysoft1'];
    end
    if wave_level==3
        cd1=detcoef(c,l,1);
        cd2=detcoef(c,l,2);
        cd3=detcoef(c,l,3);
        if cd3 > thr
        disp(cd3);
        end
        %% 进行threshold阈值处理
        ysoft1=S2GetThreshold(cd1,option,thr);
        ysoft2=S2GetThreshold(cd2,option,thr);
        ysoft3=S2GetThreshold(cd3,option,thr);
        c1 = [ca;ysoft3';ysoft2';ysoft1'];%用于高频信号
%         c1 = [ca,ysoft3,ysoft2,ysoft1];%用于低频信号
    end
    if wave_level==4
        cd1=detcoef(c,l,1);
        cd2=detcoef(c,l,2);
        cd3=detcoef(c,l,3);
        cd4=detcoef(c,l,4);
        %% 进行threshold阈值处理
        ysoft1=S2GetThreshold(cd1,option,thr);
        ysoft2=S2GetThreshold(cd2,option,thr);
        ysoft3=S2GetThreshold(cd3,option,thr);
        ysoft4=S2GetThreshold(cd4,option,thr);
        c1 = [ca;ysoft4';ysoft3';ysoft2';ysoft1'];
    end
    if wave_level==5
        cd1=detcoef(c,l,1);
        cd2=detcoef(c,l,2);
        cd3=detcoef(c,l,3);
        cd4=detcoef(c,l,4);
        cd5=detcoef(c,l,5);
        %% 进行threshold阈值处理
        ysoft1=S2GetThreshold(cd1,option,thr);
        ysoft2=S2GetThreshold(cd2,option,thr);
        ysoft3=S2GetThreshold(cd3,option,thr);
        ysoft4=S2GetThreshold(cd4,option,thr);
        ysoft5=S2GetThreshold(cd5,option,thr);
        c1 = [ca;ysoft5';ysoft4';ysoft3';ysoft2';ysoft1'];
    end
    sig=waverec(c1,l,wave_name);
    disp(thr);
end

%% 阈值函数5 -- 改进方法s3 -- S3Threshold  
function sig = S3Threshold(data,wave_name,wave_level,select_topic,option)
     %% 用wave_name小波对含噪信号进行wave_level层分解并提取系数
     [c,l]=wavedec(data,wave_level,wave_name);
     %disp(c);
     %取第wave_level层低频近似系数
     ca=appcoef(c,l,wave_name,wave_level);
     %取各层高频细节系数
     thr=thselect(data,select_topic); % 阈值获取
     if wave_level==1
        cd1=detcoef(c,l,1);
        %% 进行threshold阈值处理
        ysoft1=S3GetThreshold(cd1,option,thr);
        c1 = [ca;ysoft1'];
    end
    if wave_level==2
        cd1=detcoef(c,l,1);
        cd2=detcoef(c,l,2);
        %% 进行threshold阈值处理
        ysoft1=S3GetThreshold(cd1,option,thr);
        ysoft2=S3GetThreshold(cd2,option,thr);
        c1 = [ca;ysoft2';ysoft1'];
    end
    if wave_level==3
        cd1=detcoef(c,l,1);
        cd2=detcoef(c,l,2);
        cd3=detcoef(c,l,3);
        if cd3 > thr
        disp(cd3);
        end
        %% 进行threshold阈值处理
        ysoft1=S3GetThreshold(cd1,option,thr);
        ysoft2=S3GetThreshold(cd2,option,thr);
        ysoft3=S3GetThreshold(cd3,option,thr);
        c1 = [ca;ysoft3';ysoft2';ysoft1'];%用于高频信号
%         c1 = [ca,ysoft3,ysoft2,ysoft1];%用于低频信号
    end
    if wave_level==4
        cd1=detcoef(c,l,1);
        cd2=detcoef(c,l,2);
        cd3=detcoef(c,l,3);
        cd4=detcoef(c,l,4);
        %% 进行threshold阈值处理
        ysoft1=S3GetThreshold(cd1,option,thr);
        ysoft2=S3GetThreshold(cd2,option,thr);
        ysoft3=S3GetThreshold(cd3,option,thr);
        ysoft4=S3GetThreshold(cd4,option,thr);
        c1 = [ca;ysoft4';ysoft3';ysoft2';ysoft1'];
    end
    if wave_level==5
        cd1=detcoef(c,l,1);
        cd2=detcoef(c,l,2);
        cd3=detcoef(c,l,3);
        cd4=detcoef(c,l,4);
        cd5=detcoef(c,l,5);
        %% 进行threshold阈值处理
        ysoft1=S3GetThreshold(cd1,option,thr);
        ysoft2=S3GetThreshold(cd2,option,thr);
        ysoft3=S3GetThreshold(cd3,option,thr);
        ysoft4=S3GetThreshold(cd4,option,thr);
        ysoft5=S3GetThreshold(cd5,option,thr);
        c1 = [ca;ysoft5';ysoft4';ysoft3';ysoft2';ysoft1'];
    end
    sig=waverec(c1,l,wave_name);
    disp(thr);
end

%% 阈值函数6 -- 改进方法s4 -- S4Threshold
function sig = S4Threshold(data,wave_name,wave_level,option)
     %% 用wave_name小波对含噪信号进行wave_level层分解并提取系数
     [c,l]=wavedec(data,wave_level,wave_name);
     %取第wave_level层低频近似系数
     ca=appcoef(c,l,wave_name,wave_level);
     %取各层高频细节系数
     thr=getthr(data,wave_level,option); % 阈值获取
     if wave_level==1
        cd1=detcoef(c,l,1);
        %% 进行threshold阈值处理
        ysoft1=S4GetThreshold(cd1,option,thr);
        c1 = [ca;ysoft1'];
    end
    if wave_level==2
        cd1=detcoef(c,l,1);
        cd2=detcoef(c,l,2);
        %% 进行threshold阈值处理
        ysoft1=S4GetThreshold(cd1,option,thr);
        ysoft2=S4GetThreshold(cd2,option,thr);
        c1 = [ca;ysoft2';ysoft1'];
    end
    if wave_level==3
        cd1=detcoef(c,l,1);
        cd2=detcoef(c,l,2);
        cd3=detcoef(c,l,3);
        %% 进行threshold阈值处理
        ysoft1=S4GetThreshold(cd1,option,thr);
        ysoft2=S4GetThreshold(cd2,option,thr);
        ysoft3=S4GetThreshold(cd3,option,thr);
        c1 = [ca;ysoft3';ysoft2';ysoft1'];%用于高频信号
%         c1 = [ca,ysoft3,ysoft2,ysoft1];%用于低频信号
    end
    if wave_level==4
        cd1=detcoef(c,l,1);
        cd2=detcoef(c,l,2);
        cd3=detcoef(c,l,3);
        cd4=detcoef(c,l,4);
        %% 进行threshold阈值处理
        ysoft1=S4GetThreshold(cd1,option,thr);
        ysoft2=S4GetThreshold(cd2,option,thr);
        ysoft3=S4GetThreshold(cd3,option,thr);
        ysoft4=S4GetThreshold(cd4,option,thr);
        c1 = [ca;ysoft4';ysoft3';ysoft2';ysoft1'];
    end
    if wave_level==5
        cd1=detcoef(c,l,1);
        cd2=detcoef(c,l,2);
        cd3=detcoef(c,l,3);
        cd4=detcoef(c,l,4);
        cd5=detcoef(c,l,5);
        %% 进行threshold阈值处理
        ysoft1=S4GetThreshold(cd1,option,thr);
        ysoft2=S4GetThreshold(cd2,option,thr);
        ysoft3=S4GetThreshold(cd3,option,thr);
        ysoft4=S4GetThreshold(cd4,option,thr);
        ysoft5=S4GetThreshold(cd5,option,thr);
        c1 = [ca;ysoft5';ysoft4';ysoft3';ysoft2';ysoft1'];
    end
    sig=waverec(c1,l,wave_name);
    disp(thr);
end

%% 阈值函数7 -- 改进方法s5 -- S5Threshold  
function sig = S5Threshold(data,wave_name,wave_level,select_topic,option)
     %% 用wave_name小波对含噪信号进行wave_level层分解并提取系数
     [c,l]=wavedec(data,wave_level,wave_name);
     %取第wave_level层低频近似系数
     ca=appcoef(c,l,wave_name,wave_level);
     %取各层高频细节系数
     thr=thselect(data,select_topic); % 阈值获取
     if wave_level==1
        cd1=detcoef(c,l,1);
        %% 进行threshold阈值处理
        ysoft1=S5GetThreshold(cd1,option,thr);
        c1 = [ca;ysoft1'];
    end
    if wave_level==2
        cd1=detcoef(c,l,1);
        cd2=detcoef(c,l,2);
        %% 进行threshold阈值处理
        ysoft1=S5GetThreshold(cd1,option,thr);
        ysoft2=S5GetThreshold(cd2,option,thr);
        c1 = [ca;ysoft2';ysoft1'];
    end
    if wave_level==3
        cd1=detcoef(c,l,1);
        cd2=detcoef(c,l,2);
        cd3=detcoef(c,l,3);
        %% 进行threshold阈值处理
        ysoft1=S5GetThreshold(cd1,option,thr);
        ysoft2=S5GetThreshold(cd2,option,thr);
        ysoft3=S5GetThreshold(cd3,option,thr);
        c1 = [ca;ysoft3';ysoft2';ysoft1'];%用于高频信号
%         c1 = [ca,ysoft3,ysoft2,ysoft1];%用于低频信号
    end
    if wave_level==4
        cd1=detcoef(c,l,1);
        cd2=detcoef(c,l,2);
        cd3=detcoef(c,l,3);
        cd4=detcoef(c,l,4);
        %% 进行threshold阈值处理
        ysoft1=S5GetThreshold(cd1,option,thr);
        ysoft2=S5GetThreshold(cd2,option,thr);
        ysoft3=S5GetThreshold(cd3,option,thr);
        ysoft4=S5GetThreshold(cd4,option,thr);
        c1 = [ca;ysoft4';ysoft3';ysoft2';ysoft1'];
    end
    if wave_level==5
        cd1=detcoef(c,l,1);
        cd2=detcoef(c,l,2);
        cd3=detcoef(c,l,3);
        cd4=detcoef(c,l,4);
        cd5=detcoef(c,l,5);
        %% 进行threshold阈值处理
        ysoft1=S5GetThreshold(cd1,option,thr);
        ysoft2=S5GetThreshold(cd2,option,thr);
        ysoft3=S5GetThreshold(cd3,option,thr);
        ysoft4=S5GetThreshold(cd4,option,thr);
        ysoft5=S5GetThreshold(cd5,option,thr);
        c1 = [ca;ysoft5';ysoft4';ysoft3';ysoft2';ysoft1'];
    end
    sig=waverec(c1,l,wave_name);
    disp(thr);
end

%% 阈值函数8 -- 改进方法s6 -- S6Threshold  
function sig = S6Threshold(data,wave_name,wave_level,option)
     %% 用wave_name小波对含噪信号进行wave_level层分解并提取系数
     [c,l]=wavedec(data,wave_level,wave_name);
     %取第wave_level层低频近似系数
     ca=appcoef(c,l,wave_name,wave_level);
     %取各层高频细节系数
     thr=getthr_new(data,wave_level,option); % 阈值获取
     if wave_level==1
        cd1=detcoef(c,l,1);
        %% 进行threshold阈值处理
        ysoft1=S6GetThreshold(cd1,option,thr);
        c1 = [ca;ysoft1'];
    end
    if wave_level==2
        cd1=detcoef(c,l,1);
        cd2=detcoef(c,l,2);
        %% 进行threshold阈值处理
        ysoft1=S6GetThreshold(cd1,option,thr);
        ysoft2=S6GetThreshold(cd2,option,thr);
        c1 = [ca;ysoft2';ysoft1'];
    end
    if wave_level==3
        cd1=detcoef(c,l,1);
        cd2=detcoef(c,l,2);
        cd3=detcoef(c,l,3);
        %% 进行threshold阈值处理
        ysoft1=S6GetThreshold(cd1,option,thr);
        ysoft2=S6GetThreshold(cd2,option,thr);
        ysoft3=S6GetThreshold(cd3,option,thr);
        c1 = [ca;ysoft3';ysoft2';ysoft1'];%用于高频信号
%         c1 = [ca,ysoft3,ysoft2,ysoft1];%用于低频信号
    end
    if wave_level==4
        cd1=detcoef(c,l,1);
        cd2=detcoef(c,l,2);
        cd3=detcoef(c,l,3);
        cd4=detcoef(c,l,4);
        %% 进行threshold阈值处理
        ysoft1=S6GetThreshold(cd1,option,thr);
        ysoft2=S6GetThreshold(cd2,option,thr);
        ysoft3=S6GetThreshold(cd3,option,thr);
        ysoft4=S6GetThreshold(cd4,option,thr);
        c1 = [ca;ysoft4';ysoft3';ysoft2';ysoft1'];
    end
    if wave_level==5
        cd1=detcoef(c,l,1);
        cd2=detcoef(c,l,2);
        cd3=detcoef(c,l,3);
        cd4=detcoef(c,l,4);
        cd5=detcoef(c,l,5);
        %% 进行threshold阈值处理
        ysoft1=S6GetThreshold(cd1,option,thr);
        ysoft2=S6GetThreshold(cd2,option,thr);
        ysoft3=S6GetThreshold(cd3,option,thr);
        ysoft4=S6GetThreshold(cd4,option,thr);
        ysoft5=S6GetThreshold(cd5,option,thr);
        c1 = [ca;ysoft5';ysoft4';ysoft3';ysoft2';ysoft1'];
    end
    sig=waverec(c1,l,wave_name);
    disp(thr);
end
       
% 使用小波阈值法实现滤波，具体实现在kwden.m、kwdencmp.m、kwthresh.m、kthselect.m文件中
% 原函数分别为wden,wdencmp,wthresh,thselect，
% 输入：
% data：待滤波数据
% wave_name：小波名称
% threshold：阈值函数类型
%       threshold = 'h'时，采用硬阈值
%       threshold = 's'时，采用软阈值
%       threshold = 's1'时，采用改进方法1，改进后的软阈值：
%                     参考论文：孙万麟, 王超. 基于改进的软阈值小波包网络的电力信号消噪[J]. 海军工程大学学报, 2019(4).
%       threshold = 's2'时，采用改进方法2：
%                     参考论文：刘冲, 马立修, 潘金凤,等. 联合VMD与改进小波阈值的局放信号去噪[J]. 现代电子技术, 2021, 44(21):6.
%       threshold = 's3'时，采用改进方法4：
%                     参考论文：基于改进小波阈值-CEEMDAN算法的ECG信号去噪研究
%       threshold = 's4'时，采用改进方法3：
%                     参考论文：基于改进小波阈值函数的语音增强算法研究
%                     采用该方法需要输入两个调节因子，分别是alpha和gamma，取alpha>0，0<gamma<1
%       threshold = 's5'时，采用改进方法3：
%                     参考论文：基于VMD与改进小波阈值的地震信号去噪方法研究
%                     采用该方法需要输入两个调节因子，分别是alpha和beta
%       threshold = 's6'时，采用改进方法--自创：
%                     采用该方法需要输入两个调节因子，分别是alpha和beta
% wave_level：  小波分解水平
% select_topic： 阈值选择规则，有可选类型：
%       '' — Adaptive threshold selection using the principle of Stein's Unbiased Risk Estimate (SURE).
%       'sqtwolog' — Fixed-form threshold is sqrt(2*log(length(X))).
%       'heursure' — Heuristic variant of 'rigrsure' and 'sqtwolog'.
%       'minimaxi' — Minimax thresholding.
%       'visushrink' - 通用阈值去噪方法
% options：部分阈值函数需要补充的参数设置，即某些调节因子，需要用结构体方式调用幅值，如options.a4_alpha=2。如果不需要设置，可以赋值为options=[]。具体如下：
%          -a4_alpha：改进方法4的alpha值设置
%          -a4_gamma：改进方法4的gamma值设置
%          -a5_alpha：改进方法5的alpha值设置
%          -a5_beta：改进方法5的beta值设置
%          -a6_alpha:改进方法--自创的alpha值设置
%          -a6_beta:改进方法--自创的beta值设置
% 输出：
% s：经过滤波后得到的数据

%% 进行滤波效果评估指标计算的FilterEffectEvaluation
function [SNR,MSE,NCC] = FilterEffectEvaluation(sig_origin,sig_filter,N)
    %% 计算SNR
    Psig=sum(sig_origin'*sig_origin)/N;
    Pnoi1=sum((sig_filter-sig_origin)'*(sig_filter-sig_origin))/N;
    SNR=abs(10*log10(Psig/Pnoi1));%生成的初始加噪信号的信噪比
    
    %% 计算MSE
    error = sig_filter-sig_origin;
    MSE = mse(error);
    
    %% 计算NCC
    ncc1 = sum(sig_origin.*sig_filter);
    ncc2 = sqrt(sum(sig_origin.^2).*sum(sig_filter.^2));
    NCC = ncc1/ncc2;
    

% 对滤波效果进行评估
% 目前包含的评估指标有：信噪比SNR，均方差MSE，波形相似系数NCC
% 输入：
% sig_origin：无噪声的原始数据，一维序列
% sig_filter：滤波后的数据，一维序列
% 输出：
% SNR：信噪比
% MSE：均方误差
% NCC：波形相似系数
end

%% S4改进自适应阈值获取
function thr = getthr(data,wave_level,option)
    n = length(data);
    N = option.N;
    data_all = data(1);
    for i = 2:n
        data_all = data(i).*data_all;
    end
    ug = abs(data_all)^(1/n);
    ua = sum(data)/n;
    Gamma = (1/10)*log10(ug/ua);
    Alpha = sum((data-mean(data)).^2)/n;
    g_Gamma = sqrt(abs(2/log(0.05*Gamma)));
    thr = real(Alpha*sqrt(2*log(N))/(2^(wave_level-1)*log(wave_level+1)*g_Gamma))-(option.a4_alpha/option.a4_gamma)*(340/n);
%     thr = real(Alpha*sqrt(2*log(N))/(2^(wave_level-1)*log(wave_level+1)*g_Gamma));
    disp('阈值判断')
    disp(ug);
    disp(ua);
    disp(Gamma);
    disp(Alpha);
    disp(g_Gamma);
    disp(-(option.a4_alpha/option.a4_gamma)*(120/n));
    disp('判断结束');

%     g_Gamma = sqrt(2/log(0.05*option.a4_gamma));
%     thr = option.a4_alpha*sqrt(2*log(N))/2^(wave_level-1)*log(wave_level+1)*g_Gamma;
end

%% S6自制自适应阈值获取（创新点）
function thr = getthr_new(data,wave_level,option)
    n = length(data);
    N = option.N;
    data_new = abs(data).^2;
    adv = mean(data_new);
    data_all = data_new(1);
    for i = 2:n
        data_all = data_new(i)*data_all;
    end
    ug = data_all^(1/n);
    ua = sqrt(sum(data_new.^2))/n;
    Gamma = 10*log10(ug/ua)*(log(2)*exp(2/option.a6_alpha));
    Alpha = var(data_new);
    g_Gamma = sqrt(2/log(0.05*Gamma));
    g_Beta = exp(option.a6_beta)/(4*log(2-option.a6_beta));
    thr = Alpha*sqrt(2*log(N))/(2^(wave_level-1)*log(wave_level+1))*g_Gamma*g_Beta;
%     thr = real(Alpha*sqrt(2*log(N))/(2^(wave_level-1)*log(wave_level+1)*g_Gamma));
    disp('阈值判断')
    disp(ug);
    disp(ua);
    disp(Gamma);
    disp(Alpha);
    disp(g_Gamma);
    disp(g_Beta);
    disp('判断结束');
    
end

%% s软阈值函数
function ysoft = SGetThreshold(cd,option,thr)
    for i = 1:length(cd)
        if abs(cd(i)) >= thr
            ysoft(i) = sign(cd(i))*(abs(cd(i))-thr);
        else
            ysoft(i) = 0;
        end
    end
%     disp(ysoft(1));
%     disp('s')
end

%% h硬阈值函数
function ysoft = HGetThreshold(cd,option,thr)
    for i = 1:length(cd)
        if abs(cd(i)) >= thr
            ysoft(i) = cd(i);
        else
            ysoft(i) = 0;
        end
    end
%     disp(ysoft(1));
%     disp('h')
end

%% S1改进阈值函数
function ysoft = S1GetThreshold(cd,option,thr)
    for i = 1:length(cd)
        if abs(cd(i)) >= thr
            ysoft(i) = sign(cd(i))*((abs(cd(i))).^2-thr.^2).^(1/2);
        else
            ysoft(i) = 0;
        end
    end
%     disp(ysoft(1));
%     disp('s1');
end

%% S2改进阈值函数
function ysoft = S2GetThreshold(cd,option,thr)
    for i = 1:length(cd)
        if abs(cd(i)) >= thr
            ysoft(i) = sign(cd(i))*(abs(cd(i))-2^(thr-abs(cd(i))));
        else
            ysoft(i) = 0;
        end
    end
%     disp(ysoft(1));
%     disp('s2');
end

%% S3改进阈值函数
function ysoft = S3GetThreshold(cd,option,thr)
    for i = 1:length(cd)
        if abs(cd(i)) >= thr
            ysoft(i) = sign(cd(i))*(abs(cd(i))-(2*thr)/(exp((abs(cd(i))-thr)/(thr)+1)));
        else
            ysoft(i) = 0;
        end
    end
%     disp(ysoft(1));
%     disp('s2');
end

%% S4改进阈值函数
function ysoft = S4GetThreshold(cd,option,thr)
    for i = 1:length(cd)
        if abs(cd(i)) >= thr
            ysoft(i) = sign(cd(i))*(abs(cd(i))-thr^2/(2*abs(cd(i)))*exp(2*(thr-abs(cd(i)))));
        else
            ysoft(i) = sign(cd(i))*(thr*(exp(8*abs(cd(i)))-exp(8*option.a4_p))/(2*(exp(8*thr)-exp(8*option.a4_p))));
        end
    end
%     disp(ysoft(1));
%     disp('s4');
end

%% S5改进阈值函数
function ysoft = S5GetThreshold(cd,option,thr)
    for i = 1:length(cd)
        if abs(cd(i)) >= thr
            ysoft(i) = sign(cd(i))*abs(abs(cd(i))-thr/(abs(cd(i))^option.a5_beta-abs(thr)^option.a5_beta+1)^(1/option.a5_beta)/exp((abs(cd(i))-abs(thr))^(1/option.a5_alpha)));
        else
            ysoft(i) = 0;
        end
    end
%     disp(ysoft(1));
%     disp('s5');
end

%% S6改进阈值函数
function ysoft = S6GetThreshold(cd,option,thr)
    for i = 1:length(cd)
        if abs(cd(i)) >= thr && abs(cd(i)) <= option.a6_gamma
            ysoft(i) = sign(cd(i)).*sqrt(-2*(option.a6_sigma^2)*log(sqrt(2*pi)*option.a6_sigma*(-sign(cd(i)).*cd(i)+thr+1/(sqrt(2*pi)*option.a6_sigma))))+option.a6_mu;
        elseif abs(cd(i)) <= thr
            ysoft(i) = 0;
        else
            ysoft(i) = cd(i);
        end
    end
%     disp(ysoft(1));
%     disp('s6');
end

%% 3.ICEEMDAN分解及频谱图
function [modes]=iceemdan(x,Nstd,NR,MaxIter,SNRFlag)
    desvio_x=std(x);
    x=x/desvio_x;
    [a,b]=size(x);
    temp=zeros(b,1);
     modes=zeros(b,1);
     aux=zeros(a,b);
    for i=1:NR
        white_noise{i}=randn(size(x));%creates the noise realizations
    end

    for i=1:NR
        modes_white_noise{i}=emd(white_noise{i},'display',0);%calculates the modes of white gaussian noise
    end
    % save interval modes_white_noise
    for i=1:NR %calculates the first mode
        xi=x+Nstd*modes_white_noise{i}(:,1)'/std(modes_white_noise{i}(:,1));
        [temp, o, it]=emd(xi,'MaxNumIMF',1,'SiftMaxIterations',MaxIter,'display',0);
        aux=aux+(xi-temp')/NR;% nnnnnnnnnnnnnnnnJub局部包络
    end

    modes= (x-aux)'; %saves the first mode
    medias = aux; %  r1
    k=1;
    aux=zeros(a,b);
    es_imf = min(size(emd(medias(1,:),'SiftMaxIterations',MaxIter,'display',0)));
    while es_imf>1 %calculates the rest of the modes
        for i=1:NR
            tamanio=size(modes_white_noise{i});
            if tamanio(2)>=k+1
                noise=modes_white_noise{i}(:,k+1);
                if SNRFlag == 2
                    noise=noise/std(noise); %adjust the std of the noise
                end
                noise=Nstd*noise;
                try
                    [temp,o,it]=emd(medias(1,:)+std(medias(1,:))*noise','MaxNumIMF',1,'SiftMaxIterations',MaxIter,'display',0);
                catch    
                    temp=emd(medias(1,:)+std(medias(1,:))*noise','MaxNumIMF',1,'SiftMaxIterations',MaxIter,'display',0);
                end
            else
                try
                    [temp, o, it]=emd(medias(1,:),'MaxNumIMF',1,'SiftMaxIterations',MaxIter,'display',0);
                catch
                    temp=emd(medias(1,:),'MaxNumIMF',1,'SiftMaxIterations',MaxIter,'display',0);
                end
            end
            aux=aux+(medias(1,:)+std(medias(1,:))*noise'-temp')/NR;% r2 r3 r...
        end
        modes=[modes (medias(1,:)-aux)'];
        medias = aux;
        aux=zeros(size(x));
        k=k+1;
        es_imf = min(size(emd(medias(1,:),'SiftMaxIterations',MaxIter,'display',0)));
    end
    modes = [modes (medias(1,:))'];
    modes=modes*desvio_x;
end
% 画信号ICEEMDAN分解与各IMF分量频谱对照图
% 输入：
% y为待分解信号
% FsOrT为采样频率或采样时间向量，如果为采样频率，该变量输入单个值；如果为时间向量，该变量为与y相同长度的一维向量
% Nstd为附加噪声标准差与Y标准差之比
% NE为对信号的平均次数
% MaxIter：最大迭代次数
% 输出：
% imf为经ICEEMDAN分解后的各imf分量值
%% 主成分分析求方差贡献率函数
function per = getper(modes,ecg,a)
    modes = modes';
    ecg = ecg';
    for i = 1:1:a
        varm(i) = var(modes(i,:));%求各变量方差
        vare = var(ecg);%求原变量方差
        per(i) = varm(i)/vare;%输出方差贡献率
    end
end